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Spreading in Disordered Lattices with Different Nonlinearities 



o 



CO 



c 

o 

> 

o 

m 

(N 
O 
O 



Mario Mulansky and Arkady Pikovsky 

Department of Physics and Astronomy, Potsdam University, mi6 Potsdam, Germany, EU 



PACS 05 . 45 . -a - Nonlinear dynamics and chaos 

PACS 63 . 50 . -x - Vibrational states in disordered systems 

PACS 72. 15. Rn - Localization effects (Anderson or weak localization) 

Abstract. - We study the spreading of initially localized states in a nonlinear disordered lattice 
described by the nonlinear Schrodinger equation with random on-site potentials - a nonlinear 
generalization of the Anderson model of localization. We use a nonlinear diffusion equation to 
describe the subdiffusive spreading. To confirm the self-similar nature of the evolution we char- 
acterize the peak structure of the spreading states with help of Renyi entropies and in particular 
with the structural entropy. The latter is shown to remain constant over a wide range of time. 
Furthermore, we report on the dependence of the spreading exponents on the nonlinearity index 
in the generalized nonlinear Schrodinger disordered lattice, and show that these quantities are 
in accordance with previous theoretical estimates, based on assumptions of weak and very weak 
chaoticity of the dynamics. 



In disordered ID lattices, all eigenmodes are exponen- 
tially localized due to Anderson localization pQ. These 
models first appeared in the area of disordered electronic 
systems [2, 3 , but they are also applicable to a wide vari- 
ety of phenomena in a general context of waves (optical, 
acoustical, etc.) in disordered media [1HB]. Localization 
effectively prevents spreading of energy in such situations. 

By considering waves of large amplitudes, one faces non- 
linearity and naturally encounters the question whether 
the nonlinearity destroys the localization or not. Al- 
though this question has already been addressed numer- 
ically [7HT2]. experimentally in BECs [13TH5] and optical 
waveguides [TSKTT] as well as mathematically [TS], a full 
understanding is still elusive. It is easier to understand 
how nonlinearity destroys localization leading to thermal- 
ization [19] and self-transparency |20j in short random lat- 
tices, than to analyze asymptotic regimes at large times in 
long lattices. For the latter setups a similarity between the 
quantum kicked rotor and a ID Anderson model [21 23 
has provided an alternative realization of the effects of 
nonlinearity. 

In this paper, we study structural properties of the 
spreading field in nonlinear disordered lattices, focusing 
on their dependence on the nonlinearity index. Indeed, 
initial studies of the spreading of perturbations [7ll9"lH0] 
have been almost exclusively restricted to the behavior of 
the second moment of the distribution and of the partic- 
ipation number. These quantities, however, do not allow 



one to distinguish between all possible scenaria. In partic- 
ular, the second moment of the distribution can grow due 
to a uniform spreading of the field, but also when localized 
packages move in opposite directions. Additionally, both 
these processes may coexist with some bursts that do not 
spread at all. In order to resolve these structural features 
in a statistical way, we apply for the first time a charac- 
terization of the spreading fields in nonlinear lattices with 
generalized Renyi entropies. For guidance, we compare 
the spreading properties with that of the nonlinear diffu- 
sion equation and study the relation between the effective 
diffusion index with the nonlinearity index of the original 
model. 

Our basic model is described by the following gener- 
alization of the Discrete Anderson Nonlinear Schrodinger 
Equation (gDANSE): 



(1) 



Here n — 1, . . . , N is the lattice site index and V n is the 
uncorrelated random potential, chosen uniformly from the 
intervall [—W/2,W/2]. The coefficient f3 is proportional 
to the nonlinear strength (hereafter we assume a normal- 
ization J2 n \4>n 2 — !)• I n this work, we consider only 
the case /3 = 1 and W = 4. The parameter a, which 
we call nonlinearity index, is a novelity compared to the 
standard DANSE model with a = 1 \9§W\. Without non- 
linearity (3 — 0, eq. (fTJ) is a standard Anderson model 
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describing a disordered lattice. The Hamilton operator 
can then be diagonalized, leading to a system of eigen- 
functions $k,n with energy eigenvalues e^. An arbitrary 
wave function ip can be decomposed into these eigenfunc- 
tions rf> n = Zk C k®k,n, with C k (t) = t7 fc (0)e- fc **. With 
nonvanishing nonlinearity, this decomposition into eigen- 
functions of the linear part of the Hamilton operator is still 
possible, but now the coefficients C k are coupled through 
a nonlinearity: 



a reasonable framework for an average spreading behavior 
of localized states in the DANSE model, to be compared 
with numerical findings. 

Asymptotically, for a > 0, the spreading in © is de- 
scribed by the self-similar solution (see inset in fig. [TJ 
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The expression for the nonlinear term J\f(C) is rather cum- 
bersome (see [11] for details). In the case of integer in- 
dices a = 1, 2, 3, . . . the nonlinear coupling M(C) can be 
simplyfied by introducing an overlap matrix V of 2a + 2 
eigenf unctions. Using this, eq. ([2]) reads as: 
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For the gDANSE model ([TJ one is interested in the 
spreading of initially bounded perturbations (i.e. at time 
t = only several sites of the lattice are excited tp n ^ 
while ip n = outside). Because linear modes are ex- 
ponentially localized, these initial states are equivalent to 
an initial seeding of a finite number of modes. A qualita- 
tive picture of the spreading, based on eq. j2j , looks as fol- 
lows: the nonlinear coupling leads to chaotic dynamics of 
excited modes, as result a chaotic force acts on nonexcited 
modes and leads to their growth, and so on. Although 
several attempts to explain the detailed mechanisms of 
spreading have been made based on nonlinear eigenmode 
interactions (see, e.g., [TU] and [55]) a convincing, detailed 
description could not be found, yet. This relates to a still 
missing general understanding of statistical properties of 
weak chaos in high-dimensional Hamiltonian systems (cf. 
concept of fast Arnold diffusion developed by Chirikov and 
Vecheslavov [24]). Nevertheless, simplifying assumptions 
allowed to develop a phenomenological picture of a slow, 
subdiffusive spreading [9, 10, 22]. One of the aims of our 
work is to test these pictures by numerics. 

As the spreading is induced by the nonlinear term in 
eq. ([TJ, we suggest as a phenomelogical description the 
nonlinear diffusion equation for the probability density 
P = M 2 : 



dt dx 



ox 



(3) 



where the effective diffusion coefficient obeys a power law 
dependence ~ p a . Note, that at the moment we do not see 
a way to derive ([3]) from the DNLS {TJ directly, as for this a 
detailed theory of the microscopic chaos is needed. Neither 
do we claim that a solution of ([TJ satisfies ([3]). Our hope, 
however, is that the nonlinear diffusion equation provides 



here A is a constant fixed by the normalization condition 
J pdx — 1. The position of the edge of the spreading field 
xq has the following time dependence: 
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In order to characterize the spreading in gDANSE quan- 
titatively and to compare it with the nonlinear diffu- 
sion model, one interpretes p n — \ip n \ 2 as probability 
at site n and uses typically the mean squared deviation 
(An) 2 = ((n — (n)) 2 ) and the so-called participation num- 
ber P^ 1 = J2 n Pn- Here we suggest to also use Renyi- 
Entropies 25 as a new characterization tool: 
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Obviously, I q ^i = S = — p n lnp n is the usuual Shan- 
non entropy, while the participation number is P = eJ 2 . 

For the self-similar evolution governed by the nonlinear 
diffusion equation (0J, the variance (Ax) 2 = J x 2 p(x)dx 
as well as the Renyi entropies I q — In J p(x) q dx can 
be evaluated analytically: 
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where B(x, y) is the Beta function. In this self-similar sit- 
uation all the entropies grow with the same rate. Corre- 
spondingly, the asymptotic growth indices of the entropies 
and of the mean square displacement defined as 
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have the same value v q — ^ var = rpr^. 

The main goal of introducing Renyi entropies as charac- 
terization of the spreading is to control the peak structure 
of the field, in a similar way as these entropies are used in 
the multifractal formalism. Indeed, the parameter q deter- 
mines the sensitivity of I q on the peaks of the distribution 
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Fig. 1: (color online) Absolute square of wavefunctions |"i/v| 2 
for different nonlinearity indices and times. From top to bot- 
tom: a = 0.5, t = 10 4 (green); a = f, t = fO 5 (red); a = 2, 
t = 10 s (blue). W = 4.0, P = 1.0, N = 1024 in all three cases. 
The curves are shifted vertically for a better visibility. The in- 
set shows the self-similar solution ((4]) of the nonlinear diffusion 
equation for increasing times t (top to bottom). 

p n . Larger values of q emphasize the high peaks while for 
small values of q the background of the distribution gov- 
erns I q (the entropy Iq characterizes the support of the 
distribution). Therefore, if there are large peaks that do 
not spread (but, e.g., just drift), then the Renyi entropies 
with large q will not grow. Following the evolution of these 
entropies, we can visualize changes in the peak structure 
of the distribution. 

Fig. [1] shows exemplary averaged wave functions for 
three different values of a = 1,2,3 at times t = 10 4 , 10 5 
and 10 8 . One clearly still sees the peaked plateau even 
though these wave functions were already averaged over 
time windows and disorder realizations. 

In fig. [2] we show the evolution of Renyi entropies for the 
"standard" nonlinearity index a = 1. We stress that for 
these calculations no averaging of the wave function was 
performed, instead instantaneous entropies have been av- 
eraged over time (over time intervals between succesive 
markers in the plot) and realizations of disorder. All 
entropies with q > 0.5 show almost the same growth 
rate for large times. Entropies with very small indices 
q = 0.1, 0.25 grow slightly slower, but this is not really 
relevant: these entropies effectively measure the support 
of the distribution and are dominated by highly fluctu- 
ating exponentially decaying tails of localized eigenmodes 

Another way to characterize the peak structure of the 
distributions is to look at differences between Renyi en- 
tropies. The mostly suitable choice appears to be the 
structural entropy S s t T introduced together with the lo- 
calization entropy 5i oc in [26] : 

S atl = I 1 -I 2 = S-hxP (10) 
Si oc = h= In P. (11) 
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Fig. 2: (color online) Evolution of the entropies I q (calculated 
using logarithm with base 10) in the gDANSE {TJ with a — 
f3 — 1. The solid line has slope 0.162, it is drawn as a linear fit 
for the growth of the Shannon entropy 7i = S. 

From this definition, we can say that the entropy is 
build from a localization part and a structural part 
S = Sioc + Sstr- To see the meaning of the structural en- 
tropy, we calculate it for a uniform distribution of length 
L. In this case we find S = InP = lni and S str = 0. The 
(always positive) values of the structural entropy measure 
the relative nonuniformity of a distribution. For the self- 
similar solution (21 of the nonlinear diffusion equation the 
structural entropy, according to (|8|), is constant. 

Fig. [3] shows the time dependence of the structural en- 
tropy for different values of the nonlinearity index a ob- 
tained from numerical integration of ([1]) for an initially 
localized wavefunction (details on the numerical integra- 
tion scheme follow later in this text). One can see that 
S'str remained rather constant over time while the localiza- 
tion entropy increased as InP ~ vi lni (fig. 2]). However, 
a = 1/2 and a = 1/4 showed a stronger increase of S s ti 
than the other values of a. But compared to Sioc, which 
exhibits a clear, straight growth over time, the increase 
of the structural entropy is rather small. These findings 
demonstrate that the peak structure of the spreading wave 
function remains constant with time, while the entropies 
grow as expected from the derealization effect induced by 
the nonlinearity. The wave packet spreads uniformly and 
is not dominated by a few peaks. This supports validity 
of the nonlinear diffusion equation as a suitable model for 
the spreading in gDANSE. 

Given the nonlinear diffusion equation as a suitable 
phenomenological description of the long-time behavior 
of the gDANSE model, one still has to find the relation 
between the nonlinear term in the gDANSE ([T]) and the 
one in the diffusion equation (j3|). More precisely, the 
relation of a and a is unknown. In the literature, three 
approaches have been discussed: 

A) a = 2a which will be called strong stochasticity here 
and the resulting spreading exponent va — 0.5/ (1+ct) 
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Fig. 3: (color online) Structural entropy Sstr vs. time for 
different values of a (see label). Parameters values are 
W = 4, 1.0, N = 1024. 

was derived by Flach et al. |10j under the assumption 
of completely random phases of the eigenmode am- 
plitudes Cfc. 

B) a = 3a (weak stochasticity) gives the spreading expo- 

nent as vb — 1/(2 + 3a), this law corresponds to early 
results of Shepelyansky [22] for the standard DANSE 
model and for the quantum kicked rotor model with 
nonlinearity. 

C) a — 4a (very weak stochasticity) leads to the spread- 

ing exponent vc = 0.5/(1 + 2a), this result was ob- 
tained by Flach et al. |10j by applying some argu- 
ments on reduced chaoticity of the excited modes 
compared to Shepelyansky's model. 

Note, that in the "standard" case a = 1 these models yield 
va = 1/4, vb = 1/5 and vc = 1/6, respectively. 

To find the spreading exponents we have performed ex- 
tensive numerical simulations for different nonlinearity in- 
dices a = 1/32, 1/16, 1/8, 1/4, 1/2, 1, 2, 3. We took 10 
disorder realizations with disorder strenght W = 4 and 
lattice size N — 1024 and initialized them with a single 
excitation at one lattice site 4> na = 1. It was found re- 
cently that the energy of the state, which is a conserved 
quantity, is crucial for the spreading behavior [TO] . There- 
fore, we ensured the energy of the states to be \E\ < 1 
for all of the initial conditions by artificially setting the 
potential value to zero at the starting point V no = 0. 
Hence, we are always in the center of the energy band 
where no breather should interfere with the spreading be- 
havior. Then we ran the numerical time evolution based 
on an operator splitting and the Crank-Nicolson scheme 
for the linear part. For the time discretization we used 
a step size of At = 0.1. /3 was set to 1.0 throughout all 
simulations. This integration method is unitary and hence 
preserves the total probability within the computer accu- 
racy 10 -16 and the energy was fluctuating less than 1% 
during the simulations. These simulations were done for 
each disorder realization and we computed P, (An) 2 and 




Fig. 4: (color online) Time evolution of the localization en- 
tropy Stoc (top panel) and the the second moment (An) 2 (bot- 
tom panel) in the gDANSE model for different values of the 
nonlinearity index a (see labels in fig. [3J. Other parameters 
were W = 4, ft — 1.0, N = 1024. The lines are numerical 
fits (at the final stage of the time evolution) S\ oc ~ v^lnt and 
(An) 2 ~ t 2v ^i respectively. Fitting results are plotted in fig. [5] 

Sstr for times between t = 10 4 to t = 10 8 (10 7 for a < 1) 
and averaged over (exponentially growing) time windows. 
Finally, we fitted (An) 2 ~ t 2v ™ and P ~ t" 2 for each 
realization separately and then averaged the results over 
disorder realizations. This was repeated for each value of 
the nonlinearity index a. 

In fig. 2] the results of these simulations are shown. For 
all values of a subdiffusive spreading has been found, al- 
lowing us to obtain v vax and V2 for different values of a. In 
fig. [5] the numerical results for the spreading exponent are 
compared with the exponents derived from the nonlinear 
diffusion equation for the different assumptions A), B) and 
C). The fits for i/ var (points) and V2 (triangles) gave very 
similar results and are both close to the theoretical esti- 
mates vb and vc- The values va are clearly larger than 
the numerical results for all values of a. Note, that the 
spreading for a = 2, 3 is quite slow and so the numerical 
fits can not give very relieable results for these parameter 
values. 

As a approaches zero, all estimates A)-C) converge to 
one va,b,c(& 0) — > 1, as the nonlinear diffusion equa- 
tion becomes the usual linear one in this limit. For the 
original model, in contrast, the nonlinear term changes to 
a linear one when a — > and no spreading at all should be 
observed as the gDANSE becomes a linear equation with 
the Anderson localization property. The riddle is resolved 
by noting the role of the constant D(a) in the nonlinear 
diffusion equation ([3]). According to the Anderson local- 
ization picture, we have to set D(0) = 0, what assures no 
spreading in the linear case. This picture corresponds to 
the recent results by Veksler et al. 12 . They found de- 
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Fig. 5: Exponents of the spreading law in dependence of the 
nonlinearity index a. Circles are the values obtained from the 
numerical fits of (An) 2 ~ £ 2l/ ™r alm ^ ne triangles come from the 
fits of Sioc ~ U2\nt (compare fig. [4| . The values are slightly 
shifted horizontally for a better distinguishability. The solid 
lines are the exponents given by ([9]) for the three stochastic- 
ity assumptions: A) strong stochasticity, B) weak stochasticity, 
C) very weak stochasticity. 



creasing spreading exponents for a — > when investigating 
the short time behavior up to t — 10 [12]. A closer look 
on fig.@]also reveals that for very small a = 1/16, 1/32 the 
spreading seems to speed up between t = 10 4 and t = 10 5 . 
The graph in the right panel of fig. |6] clearifies this as one 
sees that for a = 1/16 the spreading is delayed roughly 
up to 10 4 . Furthermore, this plot examplarily shows the 
independence of our results on the time discretization At. 
Additionally, we have plotted the values of the second mo- 
ment at time t — 10 (left panel in fig. [6]) for the different 
nonlinearity indices a (note non-logarithmic scaling of a 
for a better comparability with |12j ) and the decreasing 
of (An) 2 supports the conclusion that D — >• for small 
a. One clearly sees the maximum at a = 0.25, which 
corresponds to the maximum of the initial spreading ex- 
ponent found in |12) . Our hypothesis is that the behavior 
of D{a) could be estimated via a calculation of Lyapunov 
exponents of chaos, to be reported elsewhere. 

The main conclusion of this paper is that the spreading 
of initially localized states in nonlinear disordered lattices 
can be phenomenologically well described by self similar 
solutions of the nonlinear diffusion equation. Its validity 
is supported by the finding that different Renyi entropies 
grow with the same exponent. In particular, the struc- 
tural entropy S s tv was used to measure the relation of the 
peaks and the background field in the spreading states. 
This quantity was found to remain rather constant during 
the spreading in most of the cases, supporting thus the 
self-similarity in average of the, however strongly fluctu- 
ating and highly peaked, wave function. We have shown 
the nonlinear diffusion equation to be applicable also in 
the linear limit of vanishing nonlinearity index if one as- 
sumes that the diffusion coefficient vanishes in this limit 
as well. We have found numerically that models of weak 
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Fig. 6: Left Panel: Second moment (An) 2 at time t — 10 4 vs. 
nonlinearity index a. Note the non-logarithmic scaling of the 
a-axis for a better comparability with [12]. Right Panel: Initial 
time evolution of the second moment (An) 2 for a = 1/16 and 
different values of time discretization At. 

and very weak stochasticity give good approximations for 
the spreading exponent u, but based on our numerics we 
cannot discriminate them (cf. |27|). Here, additional stud- 
ies of microscopic statistical properties of the underlying 
chaos are needed. Quite recently, it was suggested that a 
crossover between strong and very weak chaos may occur 
in the DANSE model [27] . In our simulations, we could 
not identify such a crossover. In a future work a much 
larger range of system parameters should be explored in a 
search for such an effect. 

We acknowledge useful discussions with K. Ahnert, 
M. Abel, D. Shepelyansky, S. Flach, and S. Fishman. 
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